%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 66.39 Procesamiento de Señales II
% Ledesma Francisco Hernan
% 26 de Octubre de 2008
% Version 1
% Filtro LMS (Least Mean Squared)
% TP2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inicilizacion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;

f=60;                           % Frecuencia de la señal de interferencia en Hz
fM=8192;                        % Frecuencia de muestreo
fase=60;                        % Fase de la señal de interferencia en grados 
Ai=0.3;                           % Amplitud de la señal de interferencia

%mu=0.4;                        % Paso del Filtro LMS (0<u<2/power_input)
%M=4;                           % Largo del Filtro LMS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Entrada de datos de usuario
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=input('Ingrese el largo del filtro LMS: ');
mu=input('Ingrese el paso del filtro LMS: ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generacion de las mediciones
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('C:\Documents and Settings\Hernan\Mis documentos\MATLAB\ECG_TP\señal_ECG1.mat');    
sn=ECG_signal1';                                    % Muestra de la señal pura

t=0:1:length(ECG_signal1)-1;
t=t/fM;                                             % Tiempo

fase=60/90*pi;                              % Fase de la inteferencia en radianes
w=2*pi*f;                                   % Frecuencia angular de la interferencia

in=Ai*sin(w*t+fase)';                       % Interferencia

d=sn+in;                                    % Señal afectada por la interfencia

Aj=0.1*Ai;
jn=Aj*sin(w*t)';                            % Señal de referencia

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filtro LMS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inicializacion
W0=zeros(M,1);                               
Wn=W0;                              
y=zeros(length(sn),1);
e=zeros(length(sn),1);

jn=[zeros(M-1,1); jn; zeros(M-1,1)];        % Completo con ceros la señal de referencia

% Recursion
tic;
for n=1:length(ECG_signal1)-1
    un=jn(n:n+M-1);                         % Señal de refencia en el filtro
    
    y(n)=Wn'*un;                            % Estimacion de la interfencia
    e(n)=d(n)-y(n);                         % Estimacion de la señal ECG
    
    Wn=Wn+mu*un*e(n);                       % Nuevos taps del filtro
end   
toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Graficos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mensaje1=int2str(M);
mensaje1=['Largo del filtro LMS: M=' mensaje1];

mensaje2=num2str(mu);
mensaje2=['Paso del filtro LMS: \mu=' mensaje2];

% Estimacion de la interferencia
figure(1)                                   
set(gcf,'Name','Interferencia Estimada')

plot(t(length(sn)-1000:length(sn)),y(length(sn)-1000:length(sn)))
axis([t(length(sn)-1000) t(length(sn)) -1.2 1.2])
xlabel('t[seg.]')
title('Grafico de la interferencia estimada')
text(93,1,mensaje1)
text(93,0.5,mensaje2)

% Estimacion de la señal ECG
figure(2)                                   
set(gcf,'Name','Señal ECG Estimada')

plot(t(length(sn)-10000:length(sn)),e(length(sn)-10000:length(sn)))
axis([t(length(sn)-10000) t(length(sn)) -1.5 2])
xlabel('t[seg.]')
title('Grafico de la señal ECG estimada')
text(92.6,1.5,mensaje1)
text(92.6,1.2,mensaje2)

% Errores instantaneo en la estimacion de la señal ECG
figure(3)                                   
set(gcf,'Name','Error instantaneo en la estimacion de la señal ECG')

plot(t,abs(e-sn))
axis([0 t(length(sn)) -0.01 0.25])
xlabel('t[seg.]')
title('Grafico del modulo del error instantaneo en la estimacion de la señal ECG')
text(50,0.6,mensaje1)
text(50,0.5,mensaje2)

% Señal pura real
figure(4)                                   
set(gcf,'Name','Señal pura real')

plot(t(length(sn)-10000:length(sn)),sn(length(sn)-10000:length(sn)))
axis([t(length(sn)-10000) t(length(sn)) -1.5 2])
xlabel('t[seg.]')
title('Grafico de la señal pura real')

% Señal ECG con interferencia
figure(5)                                   
set(gcf,'Name','Señal ECG con interferencia')

plot(t(length(sn)-10000:length(sn)),d(length(sn)-10000:length(sn)))
axis([t(length(sn)-10000) t(length(sn)) -1.5 2])
xlabel('t[seg.]')
title('Grafico de la señal ECG con interferencia sinusoidal')
